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NOMENCLATURE 


C Tortuosity factor. 

H Consistency index in Herschel-Bulkley model, NsVm. 

2 

k Absolute permeability of the formation, m . 

k Absolute permeability of the formation in the 

2 

presence of shear thinning behaviour of oil, m . 
k Absolute permeability of the formation in the 

n 

2 

presence of Newtonian behaviour of oil, m . 
k^^ ^ Relative permeability of oil and water respectively. 

K Thermal conductivity, W/m°C. 

L Length of the domain, m. 

n Power law index. 

P Capillary pressure between oil and water phases, 

C 

N/m^. 

2 

P^ ^ Phase pressures of oil and water respectively, N/m . 

P Absolute pressures on injection and exit planes, 

2 

N/m . 

2 

Absolute pressure at oil-water interface, N/m . 
ppv Percentage pore volume of oil recovery. 

S Saturation of oil and water in the porous region. 

O, W * w 

t, At Time and time step, s. 

T Temperature, ^C. 

T^ Temperature of rock formation, 
u Darcy velocity, m/s, 

U Composite velocity in the energy equation, m/s. 

U Velocity of the water front, m/s. 

X, y Cartesian coordinates, m. 

Ax, Ay Grid size along x and y axis repectively, m. 

Greek Symbols 
P Expansivity, ^C"^. 

? Compressibility, Pa‘\ 

p Dynamic viscosity of the fluid phase, N m/s. 

G Porosity. 

Superscript 

n Current time step. 
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Subscript 

i,J Node indices along x and y axis respectively, 
n Special case of a Newtonian fluid. 
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ABSTRACT 


A numerical study of multi-phase flow in an oil reservoir is 
presented in this work. Two phase flow of oil and water is 
mathematically modelled and the effect of surfactants in enhancing 
oil recovery is examined. Surfactants are surface active 
chemicals that reduce considerably the capillary forces in the 
fluid. Since the bulk of resistance to flow in unsaturated porous 
media comes from capillary forces, the use of surfactants can be 
expected to be quite effective in improving oil recovery. The 
governing equations are set up in terms of oil and water 
pressures. The pressure equations are seen to be coupled and 
non-linear. A fully implicit finite difference 
scheme has been used to solve the governing partial differential 
equations. Non-linearity is accounted for using Picard iterations 
within each time step. Both Newtonian as well as shear thinning 
behaviour of oil are considered in the analysis. 

In the present work a domain decomposition technique has been 
developed for solving the governing equations. Domain 
decomposition is a useful method for solving problems over very 
large domains or regions having a complex geometry. An attractive 
feature of this technique is the possibility of parallelizing the 
computer code suitable for computers with a parallel architecture. 
In the present work a completely parallel code has been developed 
for the oil recovery problem using the domain decomposition 
technique. Results do show a great potential for this technique 
while solving oil recovery problems over a field scale. 
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CHAPTER 1: INTROWCTION 


It is now widely recognized that there could be a critical world 
oil shortage in near future. In this situation a possibility of 
significantly increasing the recovery of oil from existing 
reservoirs assumes special importance. It is of course known that 
that a large proportion of the oil still remains unrecovered. 
This is so because until very recently there was no real economic 
incentive to develop enhanced oil recovery (EOR) technologies. 
Things have however changed fast and we are now beginning to see a 
sense of urgency in developing these techniques. Particularly in 
our country where oil production is negligible and EOR techniques 
are still a vague concept, much needs to be done in this 
direction. It is the objective of this work to consider the EOR 
techniques from numerical point of view and develop faster and 
more reliable numerical codes for reservoir simulation. 

Before the advent of enhanced oil recovery (EOR) techniques, 
conventional methods were used for oil recovery. These 
conventional methods fall in the following two categories: 

1) Primary recovery: In this case the pressure in the reservoir is 
high enough to force the fluid out without any pvimping effort. 
This accounts for 15%-20% of the extraction. 

2) Secondary recovery: Here water, air or steam Injection drives 
are used to maintain the reservoir pressure. Some of the oil is 
physically displaced towards the well by the high pressure 
injection. Oil extraction of around 50% is achieved at this stage. 

A large fraction of the world’s oil reservoirs remains 
unrecoverable by the conventional recovery methods. The three 
major reasons for this drawback are the following. Owing to the 
tortous nature of the microscopic pores in which oil is trapped 
only a portion of oil reserve can be contacted by the displacing 
fluid. Secondly, an Irreducable limit of oil saturation exists in 
the porous region and not all of the oil can be displaced from the 
pores by the high pressure displacing fluid. Thirdly, certain oils 
are too viscous to move at sufficient rates to support an 
economically feasible operation. Keeping these facts in mind, many 
advanced techniques collectively known as enhanced oil recovery 
(EOR) techniques, have been developed over the past two decades. 
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Field experiments have confirmed that these are more efficient 
than the conventional techniques when the oil saturation drops 
below a certain value. There are several EOR techniques in current 
use. All of them attempt to overcome one or more of the three 
factors mentioned above. For example, thermal recovery methods 
(.Boberg (1988)) aim at reducing the oil viscosity by raising its 
temperature thus improving the displacement efficiency of oil. 
Polymer flooding ( Boberg (1988), Shah (1977)) is used to increase 
the volume of the reservoir contacted by the displacing fluid. 

Lately, surfactant flooding has become popular among the EOR 
techniques because of its high efficiency and relative ease with 
which it can implemented in practice. Surfactants are surface 
active chemicals that reduce the interfacial tension between oil 
and water considerably. A surfactant solution in water injected at 
an appropriate pressure mobilizes the trapped oil to form a 
flowing oil bank by overcoming the capillary forces that are 
predominant in unsaturated porous regions. However, surfactant 
flooding is not devoid of problems. To start with, surfactants are 
very high cost chemicals that are not found in abundance. Further, 
a significant amount of laboratory work is needed to tailor the 
surfactant system for each reservoir. 

In the present work we numerically study the effect of 
surfactants on oil recovery. Oil recovery with surfactants 
dissolved in water is compared with oil recovery without 
surfactants. An appreciable increase in oil recovery in the 
presence of surfactant is seen. An ideal case is also considered 
in which the capillary resistance is zero, corresponding to a 100% 
effective surfactant. Both Newtonian and non-Newtonian behaviour 
of oil are considered throughout the analysis. For the 
non-Newtonian behaviour, oil is modelled as a power law fluid 
exhibiting shear thinning behaviour. The Herschel-Bulkley model 
(Aziz (1986)) is used to modify accordingly the governing Darcy’s 
law. Oil recovery is seen to increase further when shear thinning 
behaviour is included in the model. 

Reservoir simulation is one of the most challenging problems 
in engineering (Ewing (1983), Jim (1983)). Owing to variety of 
contributing factors like the inherent complexity of the domain, 
complex nature of the constitutive relationships and the 
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instabilty of the oil-water interface the governing partial 
differential equations are highly non-linear and coupled. While 
analytical solutions are difficult to obtain, mundane numerical 
methods exhibit instability and consequent divergence. 

In this work a pressure based formulation (Ewing (1983)) is 
used. The phase pressures are solved for using the control volume 
finite difference scheme. The scheme is fully implicit and hence 
is expected to be unconditionally stable. As mentioned before, the 
oil-water interface is inherently unstable and explicit schemes 
require excessively small time steps. The presence of mass 
transfer phenomenon such as adsorption and desorption of 
surfactants from water to the solid matrix and the vice versa Is 
neglected in the present work. The temperature of the formation Is 
taken to be a constant since the study is focussed on the effect 
of surfactants on oil recovery. 

Computational time is an important aspect to be kept in mind 
while solving reservoir problems. Because of high degree of 
non-linearity of the governing equations and large size of 
reservoirs, reservoir simulation can be computationally intensive 
even on the most powerful computers available at this time. The 
advent of parallel processing promises new generation of computers 
that can deliver an order of magnitude improvement in computing 
speeds. A parallel algorithm is accordingly needed that is 
compatible with the parallel architecture. In this context, 
domain decomposition techniques (Glowinski et al (1983)) are very 
suitable for reservoir simulation. In domain decomposition the 
physical domain is split into sub-domains. The governing 
equations are then solved Independently in each sub-domain. The 
individual solutions must be assembled carefully to get a 
converged solution of the original problem on the complete 
physical domain. Domain decomposition can rapidly make the 
simulation parallelizable since calculations in each sub-domain 
can be carried out independently on a proceessor. They can also 
handle large problems on sequential machines more efficiently. 

A large body of works on domain decomposition have appeared 
in the literature in the past few years. Many of them originate 
from the Schwarz method (Glowinski et al (1983)) for solving 
elliptic problems. This method is originally sequential in 
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nature. Lately it has been modified in view of its good 
parallelization properties (Glowinski et al (1983)). Glowinski et 
al (1983), in their pioneering work have developed a number of 
efficient parallelizable domain decomposition algorithms and these 
are currently in wide use. Perng and Street (1990) have applied 
domain decomposition technique to solve Navier-Strokes equations in 
complex domains which are otherwise difficult to solve without 
using coordinate transformation. Parralelizable domain 

decomposition codes applied to transient problems are however 
scant in literature. Sequential domain decomposition codes 
applied to steady problems are available in abundance (Perng 
(1990)). 

In the present work we have developed a parallel code to 
solve the transient oil recovery problem. Uzawa’ s algorithm 
(Glowinski (1983), Yagawa (1991)) has been used for interface 
treatment in conjunction with an implicit formulation. The 
overall algorithm has the desired parallel and temporal 
properties. Results are in excellent agreement with those of a 
full domain simulation with no observed deterioration in the 
sequential CPU time. 

The thesis is organized along the following lines: 

1) Mathematical model and formulation of the oil recovery problem. 

2) Developement of domain decomposition technique. 

3) Results of simulation of enhanced oil recovery. 

4) Scope for future work. 
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CHAPTER 2: MATHEMATICAL MODEL AND FORMULATION 


We start this section with the derivation of the equations 
governing oil and water pressures in two phase flow through a 
porous medium. 

2. 1 ASSUMPTIONS 

The following assumptions have been made: 

1) The two fluid phases namely, oil and water flow simultaneously. 

2) The fluids are immiscible i.e. there is no mass transfer 
between them. 

3) Gravity term is neglected in the momentum equation. It is 
reasonable because the time frames are short and as a result the 
gravity override effect will not show up. 

4) The porous formation is incompressible and porosity is a 
constant. 

5) Relative permeabilities of phases and capillary pressure are 
only a function of saturation. The fluid viscosities depend on 
temperature only. 

6) There is no net mass source in the flow field. 

2.2 PHYSICAL PRINCIPLES 

Following are the basic equations that govern flow through a 
homogeneous, isotropic porous region. These are used to derive 
the final pressure equations. In this derivation oil and water are 
assvuned to exhibit Newtonian behaviour. Modifications arising from 
non-Newtonian behaviour of oil are presented later, 
a) Mass balance 

Mass balance equations for the oil and water phases are the 
following. 
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oil: 


+V. u p = q 

^ o o o 


( 2 . 1 ) 


a (eS p ) 

o o 


at 


water: 


a (eS p ) 


at 


+ 7. u p = q 

^ W %f %f 


( 2 . 2 ) 


Since the formation is incompressible, e, the porosity is a 
constant. The mass source terms, q and q are taken to be zero 

O H 

in this work. p and p are the phase densities. S and S are 

o w o w 

oil and water saturation repectively. u and u are the Darcy 

O W 

velocities of oil and water respectively. We assume that no gas 
phase is present in the formation and so, 

S + S =1 (2.3) 

O 

b) Momentum equation: 

The momentum equation is taken in the form of a generalized 
Darcy’s law as follows. 


oil: 


water: 


k k 

ro 

k k 

rw 


( 7 P ) 

^ ft 


( 7 P ) 

^ W 


(2.4) 


(2.5) 


Here, k is the absolute permeability of the porous formation, k 

ro 

and k are the relative permeabilities of oil and water 

rw 

respectively, u and u are the phase viscosities and P and P 

o w o w 

are the phase pressures, 
c) Equations of state: 

The quantities ^ and P defined below are taken to be specified. 

1 

compressibility: ? * 

“ ap 


o 


T 


( 2 . 6 ) 




w 



T 


Volumetric thermal expansion coefficient: 
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ax 


(2.7) 


p = - 



p 

o 




w 



sp, 

ax 


p 

w 


d) Constitutive relations: 

Xhe formation properties given below are also taken to be 
available from laboratory or field experiments. 

Capillary pressure: 

P - P = P (S ) (2.8) 

o w c w 

Relative permeability: 

k = k (S ) (2.9) 

ro ro w 

k = k (S ) (2.10) 

rw rw w 

Here subscripts o and w stand for oil and water respectively. X 
in the definition of 3 is the local (cell -averaged) temperature. 

2.3 GOVERNING DIFFERENXIAL EQUAXION 

We incorporate the equations of state and the constitutive 
relations into the mass balance equation for oil. Using Equations 
2.1, 2.3 and 2.7, we get after simplification the following oil 
pressure equation: 


eS p 

o o 




ap 

0 

"aT 


- P 


ax 

at 


as k k p 

+ ep » 7. 7 p (2.11) 

©at ~ n ~ o 


Similarly for water pressure, we get 
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€S p 

w w 


dP 


at 


81 

at 


as k k p 

+ ep = V. 7 P (2.12) 

W Ot ~ P ~ w 


Equations 2.11 and 2.12 may be written as. 


-S 0 - 

O O 


ai 

at 


ap 

s C + 

o o ot 


as 

o 

at 


= — — r V. 

ep I ~ 


k k p 
ro o 


7 P 


(2.13) 


-S 0 

w w 


81 

at 


ap 


+ s c 

w w 


at 


as 


at 


= -i— [ 7. 
ep [ ~ 


k k p 
rw w 




(2.14) 


At this stage we introduce the capillary pressure, P , since 

C 

® ^ unique functions of P , Using Equation 2.8 

ro rw o w c 

and the relation, 

as dS ap 

w _ w c 

“at~ "dP at~ 

c 

Equations 2.13 and 2.14 reduce to the following: 
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oil: 


-S ^ 


o o 


dT 

at 


+ s ? 

o o 


SP 

at 


dS 

H 

dP~ 


c 



ap 

w 

W 



k k p 
ro I 



(2.15) 


„T 3P 

water: -S ^ + S ? 

w H at w w at 


dS 

M 

"dF 


ap 


at 


ap T 

-at^] 


1 r ^ ^ P 1 

= — I — 7. I 1 7 P 

~ L J ~ H 


(2.16) 


Equations 2.15 and 2.16 are the final form of the pressure 
equations used in the present work for numerical simulation. Once 
P^ and P^ are known other quantities can be determined from the 
constitutive relations. The energy equation in its final form is 
given as, 


+ u. 7 T = 7^ T (2. 17) 

ct (T 

T 

where U = — ^ ZpCu ,cr = e-|ESpC 1+ (1-e) (pC) 

~ <r^ 1 1 1 ’ T T I'^l 1 J ^ for»atlon 

in which Cj is the specific heat of phase i (« o, w) and K is the 
thermal conductivity of the formation. Other symbols hold the same 
meaning as mentioned before. 

It is necessary to Justify the use of pressure formulation at 
this stage. Oil recovery problem can also be formulated in terms 
of water saturation. This approach can however lead to 
difficulties in numerical computation since the oil-water 
Interface can be quite sharp. In contrast to this the pressure 
field must be continuous everywhere. This results in considerable 
improvement in the accuracy of the numerical simulation. 
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2.4 NON - NEWTONIAN BEHAVIOUR 

We look at the non-Newtonian behaviour of oil. Field experiments 
show that oil can be modelled as a power-law fluid exhibiting 
shear-thinning behaviour. Darcy’s law, as it appears in Equations 
2.3 and 2.4 may be written in its general form as 

k( y P ) k ( P ) 

u= — ( VP^- G ) (2.18) 


where the subscript i refers Individually to the oil and water 
phases. G is the thresold pressure gradient below which there is 
no flow. It has been taken to be zero in this study. We see that 
k, the absolute permeability, depends on the local pressure 
gradient. This dependence is absent when the fluids exhibit 
Newtonian behaviour. The shear-thinning behaviour of oil is 
well-described by the Herschel-Bulkley model [Aziz (1986)). In 
this model the appropriate form of k in Equation 2.18 for the oil 
phase is given as, 

k 

^ f — 7 ) — 1 I ^ P ( 2 . 19 ) 

off o I U J " o 


where k^ is the absolute permeability of the formation in the 
absence of non-Newtonian behaviour and u is given as, 

ef f 


= 

©ff 


H 

4 




(1 - n)/ 2 


( 8 C e 


k ) 

n 


( 2 . 20 ) 


In the above equation H is the consistency index of the 
Herschel-Bulkley model and n is the flow behaviour index 
equivalent to the index in a power-law fluid. C is the tortuosity 
factor of the formation. Based on the data presented by (Aziz 
(1986)), C is taken as 25/12 in the present study. Equations 2.19 
and 2.20 reduce to the Newtonian case of u = H = u and k » k 

eff o err n 

for n=l. The surfactant solution can also exhibit non-Newtonian 
behaviour. In this study the surfactant solution is however 
assumed to be Newtonian. 
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2.5 GEOMETRIC MODEL 

The geometry of the physical region used in the present study is 
shown in Figure 2.1. The flow is two dimensional and the 
rectangular domain has a nominal size of 4m x 0.4m. Thus, if the 
domain has a length of Xm and breadth of Ym then 
X = 4m and Y = 0.4m 

Further, in a Cartesian x-y plane x=0 and x=X are the inflow and 
outflow planes respectively. On these planes constant pressures 
are specified. The upper and the lower edges of the domain i.e. 
y=0 and y=0.4 are impermeable. We solve for the oil and water 
pressures separately. The difference between the two pressures is 
the capillary pressure P , and is a measure of the surface 

C 

tension. The capillary pressure prevents equalization of the 

phase pressures and hence, P^= P^- P^is non-zero. Typical values 

of the parameters appearing in Equations 2.6 and 2.7, the 

equations of state, are taken from Boberg (1988) and summarized in 

Table 2.1. The constitutive relations for an unsaturated porous 

medium namely, k = k ( S ) and P = P ( S ) are shown in Figures 
rl rl w c c *# 

2.2a and 2.2b. Since E 1, the fvmctional relationship can be 
shown in terms of one phase saturation alone. In the present 
study results have been obtained in terms of water saturation. An 
Increase in with time clearly shows oil displacement. Since the 
effect of surfactants is to reduce the capillary forces and 
mobilize the oil bank, the capillary pressure P decreases and the 

C 

relative permeabilities, k^^. Increase in its presence. This 
effect is clearly shown in Figures 2.2a and 2.2b. Curves a, b and 
c refer to the relationships in the absence of surfactant, in the 
presence of surfactant and in the presence of an ideal surfctant 
respectively. Owing to lack of availability of real field data, 
these constitutive relations show the general effect of 
surfactants. They are in tune with the relations presented in 
Barenblatt et al (1986). 

2.6 INITIAL CONDITION AND KXJNDARY CONDITIONS 

The pressure equations. Equations 2.15 and 2.16, are numerically 
solved on a rectangular domain, Figure 2.1, subject to a quiscient 
initial condition, impermeable confining walls on which 5P^/9y » 
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0, where i refers individually to the oil phase and water phase, 
and specified oil and water pressures on the inflow and outflow 
planes. This model predicts the movement of the water front from 
the inflow plane into the oil-rich porous zone. This increase in 
water saturation is a measure of the amoxint of oil displaced. Oil 

displacement is measured in terms of the percentage pore volume 

, s T 4 . 1 j total volume of oil produced 

(ppv). It is defined as, ppv = 100 x r— t-t i — 

total pore volume 

The use of Dirichlet boundary condition on the outflow plane 
is justified only if the location of this plane is far away from 
the region effected by fluid injection. In the present study oil 
recovery (ppv) is computed over a distance of 2m and the outflow 

plane is located at a distance of 4m from the inflow plane. For 

the time period of three hours considered in the present 
simulation the water moves a distance of approximately Im. 
Numerical experiments show that oil recovery is insensitive to the 
location of the outflow plane for the time period considered here. 

2.7 CASE OF IDEAL SURFACTANT 

We now dicuss the limiting case of an ideal surfactant which when 
added to the injected fluid eliminates surface tension altogether. 
Hence the relative permeabilty and capillary pressure reach their 
limiting values of unity and zero respectively as shown in Figures 

2.2a and 2.2b. A water front moves through the porous medium under 

these conditions displacing oil completely. We assume that the 
water saturation is at its maximum value ahead of the front. 
Similarly, beyond the front water saturation is minimum. Since 
the effect of surfactant is totally negated flow resistance arises 
from the phase viscosities alone. For a domain of length L the 
average oil and water velocities are equal and equal to the speed 
of the front movement. For the Newtonian case this speed is 
estimated from the principle that the front speed computed from 
the oil side be equal to that computed from the water side. 

Hence , 


k ( P - P ) 

U = — 

front fi X 


k (p.- p^) 

u ( L - x) 


( 2 . 21 ) 


Here P^, P^ and P^ are the injection plane, interface and exit 
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plane pressures respectively, x indicates the depth penetrated by 
the water-oil interface. Eliminating in Equation 2.21 the 
expression for the front speed is obtained as, 


U 


( p,- p,) 


front 


+ M ( 1 - 


X 

L 


( 2 . 22 ) 


At t = 0, X = 0 and at any time t we update x as, x = x + 

^ neH old 

U At . Here U, ^ is calculated from Equation 2.21 using the 

front o front 

old value of x. Oil recovery Is then measured in terms of ppv as, 
lOOx 

ppv = . 

When the shear-thinning behaviour of oil is taken into 
consideration, a closed form solution for U does not exist. In 
terms of P^, P^ and P^ as well as the power-law index n, the 
modified form of Darcy’s law , Equation 2.18, can be written as, 


r K T n r P - P^ Tn k rP-P n 


u 


front 


For each x Equation 2.23 is solved for P^lteratively. 
value of P the interface velocity U 

I ■' front 


(2.23) 

With this 
is calculated as. 


_ 1 _ 

n 


U 


-I pP— P •% 

•[-fc;) 


n 


front 


(2.24) 


The same procedure as given after Equation 2.20 is adopted to 
obtain X at any time t and finally the ppv. 
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Table 2.1 : Properties Used in Numerical Calculation 


Absolute Permeability, K = 132 Darcies 

Porosity, c = 0.375 

Compressibility: oil, = 0.03447 Pa"^, water, = 0.02137 Pa 

In-situ formation pressure * 1.31 MPa 

Injection pressure - 1.793 MPa 

Dynamic Viscosity (Pa-s) 


T, °C 

20 

40 

60 

80 

Oil 

4.12 

-3 

2.37 

-4 

0.61 

-4 

0. 10 

Water 

1.03x10 ^ 

6.6x10 

4.73x10 

3.67x10’ 









d 

o 


•r-i 

CD 

cts cn 



Figure 2.2b Capillary pressure as function of i 
a) Without surfactant b) With surfactant c) Ide: 



CHAPTER 3: INTRODUCTION TO DOMAIN DECOMPOSITION 


3.1 OPENING REMARKS 

Practical engineering problems Inevitably involve complex 
geometries and large domains. Computation of three dimensional 
viscous flows, high speed gas flows, turbulent flows with 
separation and enhanced oil recovery require enormous amount of 
computing time and memory space. In many case the computations 
cannot be Implemented on existing computers at all. The current 
approach towards complex problems is to improve the numerical 
algorithm on one hand and computer performance on the other hand. 
With the advent of parallel computers there exists a great 
potential in increasing computer speeds as well as memory. 
Parallel computers require parallel algorithm to exploit the 
distributed nature of the processor architecture. Domain 
decomposition techniques are the simplest of parallel algorithms 
that can convert a sequential algorithm to one suitable for 
parallel computing. Besides the ability to parallelize, other 
specific advantages of domain decomposition are given below. 

1) It can generate smaller matrices that require less memory and 
are rapidly invertible when Iterative techniques are used. 

2) Complex assemblage of components of a system can be analyzed 
systematically. 

3) Complex phenemena can be localized by assigning separate 
sub-domains to them. 

3.2 MODEL PROBLEM 

We demonstrate domain decomposition techniques first through a 
simple example. Consider a linear governing differential equation 
of the type 

Au = f in Q (3.1) 

and u = f on 90 

where f and g are known functions. Here O is the physical region 
and 9Q is its boundary as shown in Figure 3.1. A is the Laplacian 
operator. 
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Figure 3. 1 


The basic theme of domain decomposition is that a large 
domain is divided into many sub-domains that are easier to handle 
and are linked at the interfaces. These sub-domains may be 
overlapping or non-overlapping. We first consider the case of 
non-overlapping sub-domains. 



Figure 3.2 shows a division of the full domain into two 
sub-domains, 1 and 2. The interface between the two sub-domains 
is denoted C. We seek a solution of Equation 3. 1 in the original 
large domain. The problem in region 1 is solved numerically with 
an arbitrary Dlrichlet boundary condition on C. Let this solution 
be uj , u being the original dependent variable. This is followed 
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by numerical solution of u in region 2 with a boundary condition 

on C which is derived from u^ . The global solution u and its 

derivative should be continous at the interface. Hence the most 

obvious choice of the boundary condition on C for region 2 would 

be the normal derivative 3u^/9n. Here n denotes the normal at C 

1 

exterior to region 1. With this boundary condition u can be 
solved in region 2. Let this solution be u\ This marks the end 
of one complete sweep of the original domain. Thus, after one 
complete sweep of the original domain we have solved the following 
equations: 

In region 1, 


In region 2, 


Au = f 
1 

u = g on an and u = A on C 
1 ® 1 


(3.2) 


Au = f 
2 

u = g on an and 
2 




on C 
C 


(3.3) 


where A is the Dirichlet boundary condition described before and 
n^ and n^ are the normals at C exterior to regions 1 and 2 
respectively. 

The problem in region 1 is solved again with a new boundary 

condition on C. This may be the value of u^ on C or some 

2 

11 2 
combination of u and u on C. Let the solution be denoted u . 

2 1 1 

The procedure continues in this fashion where the problem In 
regions 1 and 2 are solved alternately until convergence in the 
global solution is achieved. This convergence criterion will 
satisfy the conditions of continuity of u and its derivative at 
the interface. This method described above is known as the 
alternating Dirchilet-Neumann (D-N) method because Dirchilet and 
Neumann boundary conditions are applied alternately on the 
interface. 
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Figure 3.3 

We now look at domain decomposition with overlapping 
sub-domains. The method described below is also known as the 
Schwarz-alternating method. Figure 3.3 shows two sub-domains 1 
and 2, with a slight overlap that is shaded. C^is the interface 
of region 1 and is that of region 2. We consider the solution 

of Equation 3.1 again in the full domain £1. The sequence of 

calculations is as follows. The problem in region 1 is solved 

numerically with an arbitrary boundary condition on C . Let the 

1 ^ 
solution be u^. Then the problem in region 2 is solved with a 

boundary condition on that is already obtained in region 1 i.e. 

u^ on C , since C falls in region 1. Let the solution be denoted 
2 2 

u\ Subsequent iterations are generated as follows. 

The problem in region 1 is solved, with a new boundary 

condition on C that may be the value of u^ on C or a linear 

^11 2 1 

combination of u and u on C . This is equivalent to the use of 
12 1 

under-relaxation factor for the numerical scheme. The iterations 

continue until convergence is achieved. The convergence criterion 

I I 

should ensure that after I iterations, u^ and u^ are sufficiently 
close. In a finite difference scheme if the overlap is such that 
there are a few nodes lying in it then the convergence criteria 
should be such that after I iterations, u^ and uj at these nodes 
are sufficiently close. Only then do we get a globally converged 
solution for u. 

One major advantage of Schwarz-alternating scheme is that it 
can do away with Neumann boundary condition altogether. This 
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however requires that the overlap be chosen judiciously. In a 
finite difference scheme if the overlap is only of a grid size 
then by solely using Dirichlet boundary conditions on and we 
can ensure the continuity of the derivatives at the overlap 
region. This is a major advantage because Dirichlet boxondary 
conditions are easier to handle and are computationally faster 
when iterations are involved, as in the non-linear problems. 

The two domain decomposition methods described above are 

sequential in nature. Their efficiency depends on how the 

iterations are handled at the sub-domain interfaces. Consider the 

alternating D-N method as an example. After one complete sweep of 

the domain, the use of arlthmatic mean of u^ and u'^ as the 

1 2 

boundary condition on C for u^*' can lead to a considerable 
improvement in overall convergence rate. The alternating Schwarz 
method also shows similar improvement. Thus, interfaces are the 
most crucial zones and they should be handled with care and 
intelligence. 

The ultimate utility of domain decomposition techniques lies 

in their ability to permit parallel computing. Domain 

decomposition techniques have been rediscovered primarily for this 

feature. We now describe a domain decomposition algorithm that 

permits parallelization. It Is known as Uzawa’ s algorithm. It 

employs non-overlapping sub-domains. For definitiveness the 

problem in Equation 3.1 is reconsidered. For a globally converged 

solution we require both u and its derivative to be continuous at 
the interface C (Figure 3.2). Keeping this in mind we start with 


a guess value of the derivatives, 


9u 


an 

1 


and - 


au 


an 


where n and 
1 


n^ are the normals at C exterior to regions 1 and 2 respectively. 
Let this guess be X^. Using as the boundary condition on C we 
solve the governing differential equation in regions 1 and 2. Let 


the solutions be u^ and u^. 

1 2 

u repectively on C. Then X is updated as 

^ ..2 .. 1 ^ , 1 I . 

X = X + p( y^ - y^) 


Let y^ and y^ be the values of u^ and 
■'i •'z 1 


C3.4) 


where p is a suitable constant that lies between 0 and 1. 

Using X^ as the next guess we repeat the earlier steps. Thus the 
overall algorithm has the following structure: 

1) Assume X^, where X^ is the normal derivative at the interface 



No. <^' 4 . 
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exterior to the sub-domains. 

2) Solve the governing partial differential equation in the two 
sub-domains, 1 and 2. i.e. 

In sub-domain 1 

Au = f - (3.5) 

1 5u 

u = g on 3Q and — s = X on C 



where n and n are the normals at C exterior to the sub-domains 1 
1 2 

and 2 respectively. 

3) Using the values of u^ and u^ at the interface i.e. and y^ , 

12 12 

X^ is updated as in Equation 3.4. Thus after I iterations, 

x^*^= x^ * pi yl - yj) (3.7) 

It can be seen that this scheme is completely parallel. For 
example, calculations in regions 1 and 2 can be carried out 
independently on two different processors. The processors return 
the values of y^ and y^ after each iteration. Iterations can be 
stopped when 

< c (3.8) 

where e is the desired convergence criterion. 

As mentioned before, Uzawa’s algorithm can be implemented on 
a parallel computer. We now describe briefly how this can be 
implemented on a hypercube-type parallel computer (Ganagi et ai 
(1991 )). The hypercube architecture has been adopted by a variety 
of research groups all over the world. We describe the 
implementation of the Uzawa’s algorithm on PACE, a parallel 
computer of the hyper cube- type developed by ANURAG, a defence 
research lab under the ministry of defence. Government of India. 
PACE-8 has eight concurrent processors and similarly PACE-128 has 
128 concurrent processors. Hence with PACE-8 as many as eight 
sub-domains can be handled. 

The complete parallel computer code will have two parts. The 
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first part is the computation-intensive portion that deals with 
the numerical solution of the governing equations in each 
sub-domain. These calculations are performed on individual 
processors. The second part is the front-end-processor (FEP) 
which deals with keeping track of interface values, checking for 
their convergence and downloading the initial and boundary 
conditions onto the sub-domains at the beginning of each 
iteration. At the beginning of each time step the 
front-end-processing part supplies the required initial values at 
the interface as well as the internal points of the sub-domains. 
The boundary conditions are also supplied. Then the computation 
Intensive part takes over and with the supplied values 
calculations are performed in the sub-domains on the respective 
processors in parallel. After local convergence is achieved in 
all the sub-domains, the front-end -processing part collects the 
final solutions for the sub-domains and checks the interface 
values for convergence. In this fashion iterations proceed and 
upon interface convergence the front-end-processing part moves on 
to the next time step. 

3.3 SPECIFIC ADVANTAGES OF DOMAIN DECOMPOSITION 

We list below the advantages of domain decomposition techniques in 
the solution of engineering problems. 

1) Parallelization: As discussed above, domain decomposition has 
the ability to parallelize the computer code. Hence the 
computations are faster and the problem of memory is circumvented. 

2) Complex domains: Domain decomposition techniques are much 
simpler to use in domains with complicated geometries. Using 
domain splitting the complex domain can be divided into many 
regular sub-domains. The governing equations can be solved easily 
on the regular sub-domains. 

3) Large domain simulation on sequential machines: 

Even in the absence of parallel computing domain decomposition can 
still be appreciably faster than ordinary simulations for large 
domains. Let us elaborate on this point. We consider a domain 
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decomposed problem with two sub-domains against an ordinary (full 
domain) simulation. For a finite difference code with a matrix 
size N (number of discretizedequations) the convergence rate of an 
iterative matrix varies as 1/N. Thus, if e is the residual error 
after I iterations, e is given as 

e = K 1 “'^ (3.9) 

where K and a are constants. Note that the convergence rate goes 
to zero as N goes to infinity. Hence Iterative inversion of large 
matrices is slow unless an extremely good initial guess is 
available. 

Let E be the desired level of residual error and be the 
number of Iterations required. Equation 3.9 can be written as, 

log E = log + log K (3.10) 


and so 




M/a 


Let us now split the original physical domain into two equal 
sub-domains. The order of matrix arising from finite differences 
in each sub-domain becomes N/2. The number of iterations required 
to the reach the same level of residual error , E, would now be 


I 


2 



(3.11) 


Let D be the number of Iterations required at the Interface. 
Keeping in mind that there are two sub-domains, the total number 
of Iterations for the domain decomposed problem would be 


Thus, 



(3.12) 


(3.13) 


A plot of (I /I ) versus N is shown in Figure 3.4 . 
si 
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1 

I /I 

S 1 


N 

«= N 

Figure 3.4 

Note that above a critical value of N, (I /I ) < 1, l.e. the 

s 1 

domain decomposed problem requires lesser number of iterations 
than the full domain problem. 

3.4 DOMAIN DECOMPOSITION FOR NON - LINEAR PROBLEMS 
Throughout the earlier sections we have assumed that the governing 
problem is linear (Equation 3.1). In case of a non-linear problem 
the Interface treatment must be modified. Let us take up an 
example to illustrate this point. Consider the linear heat 

conduction problem in the same domain as in Figure 3.2: 

K V^T = pc (3. 14) 

p ot 

where K, p and c are constants. 

p 

In Equation 3.14 the boundary conditions to be matched at the 

interfaces of the sub-domains would be the temperature and the 

dX 

temperature derivative, —g~ • This is both physically and 
mathematically conceivable. Consider the non-linear problem, 

V. K 7 T = pc -1^ (3.15) 

where K is taken as a function of temperature. The physically 
conceivable Neumann bovindary condition to be matched at the 
interface would be not the derivatives of T but the fluxes at the 
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interface , -K(T) — . In case of the linear problem this boils 

on 

down to the temperature derivatives alone since K is not a 
function of temperature. Thus for the non-linear problem, global 
convergence requires 


-K( tJ 



K( T 



( 3 . 16 ) 


It should be noted that during Iterations * ^2 ‘ 

inhomogeneous media with sharp change in thermal conductivity at 
the interface, heat fluxes rather than the temperature gradients 
should be matched. 
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CHAPTER 4: NUMERICAL SIMULATION OF ENHANCED OIL RECOVERY 


In the present chapter we develop the formulae applicable for 
numerical simulation of oil recovery. We first discuss the full 
domain simulation and then proceed to the formulation with domain 
decomposition. The governing differential equations and boundary 
conditions are the same as derived in Chapter 2. 

4. 1 FULL DOMAIN SIMULATION 

A central difference scheme is used for the numerical solution of 
governing differential equations. The physical domain is shown in 
Figure 4.1. The control volume is shown below, Figure 4.2. 



Figure 4.2 

Let N and M be the number of nodes in x and y directions 
respectively. Ax and Ay define a Cartesian mesh and At is the 
time step. The current level is represented by n and the future 
time level by n+1. The scheme is fully implicit with respect to 
time stepping and is expected to be unconditionally stable. 

The oil pressure equation, Equation 2.15, Is reconsidered 
here. It is given as. 
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(4.1) 

We now discretize this equation using finite differences. Since 
the formation temperature is taken to be constant in this study, 
the first term in the left hand side namely, -S ^ -- ir— » Is 

O O C/ 1 

dropped. Let 

K(VP ) K (S ) p (T) 

~ o ro H o 


Then the right hand side is 



V. 


0 7 P and is equal to 

~ o 



1 


P 


ol, j 


r P 

- p 

T n+ 1 

r P - P 

1 n+ 1 

1 ol+l,J 

ol, J 

■ •- 

1 ol,J 

1 

L 

[ Ax 

J 


Ax 


+ 


1 


P 


oJ, J 


P — p 

r oi.j-M oi,j 

1 

- a 

8 

r P - P 1 

ol,J ol,J-l 

n-i-X 

L J 


Ay 


Here 0,0,0 and 0 are evaluated as harmonic averages of their 

own s 

values at nodes that lie on either side, Figure 4.2. Hence 


6 

e 


2 0 



a 


i, J 


+ 0 
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Similarly, 


8 = harmonic average (0,8 ) 

w ® 1. l 1-1. 1 




II 

C 

CD 

harmonic average ( 6 ^ j 

, 0 ) 
i,j+i 



0 = 
s 

harmonic average ( 0j j 

, 8 ) 
l.J-1 

The left 

hand 

side 

is discretized as. 


f S ? - 

o o 

dS 

W 

.1 

r P""' - P*' , r 

1 oi.j oi,j 1 . r 

dS T r P"*^ - P" , 

W 1 1 oi.J 1 

dP 

c 

M.j 

L J [ - 

dP^ J L At J 


The full equation in discretized form may be written as. 


A P 

o Ol+l,J 


B P 

o Ol,J 


C P 

O Hl,J 


D P 

O 


E P 

o ol,J+l 


F P 

o o 




n+1 


= G 


(4.2) 

The water pressure equation. Equation 2.16, can be discretized 
similarly to yield. 


] n+l 


= G 

(4.3) 

The coefficients A - G , A - G are given in Appendix A2. 

O O W H 


The discretized boundary conditions are as follows: 
Specif led inf low plane and outflow plane pressures: 

P (l.J) = P , 1 

O ol 

P.Cl.J) = P„ 

P^(N.J) - P^^ 

P_(N,J) = P^ 


(4.4) 
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Impermeable confining walls at y = 0 and y = Y: 


P (i,l) - P (i,2) = 0 

o o 

P (1,1) - P (i,2) = 0 

w w 

P (i.M) - P (i,M-l) = 0 

O O 

P (i.M) - P (i,M-l) = 0 

H V 


(4.5) 


The pressure equations, Equations 4.2 and 4.3, are coupled 
and are solved simulataneously. When Equations 4. 2-4. 5 are 
written for each node, a set of 2NM non-linear algebraic equations 
is obtained. 

Non-linearity has been primarily handled using Picard 

iterations within each time step. In this scheme the 

coefficients, A^-G^, A^-G^, depend on P^ and P^. As a first guess 

they are evaluated using the previous time step values of P , P 

o w 

and S^. Iterations continue within each time step until 

convergence in P and P are achieved. The matrix structure 

O W 

arising from Equations 4. 2-4. 3 is shown in Figure 4.3. It is 
banded and has an overall size of 2NM x 2NM. This matrix is 
Inverted using an efficient sparse matrix inverter (Duff ( 1980 )). 


4.2 ALGORITHM 

We now outline the algorithmic steps required to obtain P^, P^ , 
and ppv as a function of time. 

1) The initial distribution of P , P and S at t = 0 are 

o w w 

prescribed. 

2) The coefficients, A -G , A -G of the pressure equations are 

O O W It 

computed using the previous time step values of P^, P^ and . 

3) The system of pressure equations is solved to obtain new values 
of P and P . 

o w 

4) From P and P , S is obtained using the constitutive 

O H W 

relations. 

5) Steps 1-4 are repeated until convergence of P^ and P^ is 
achieved within the time step. 

6) S is updated and ppv is calculated using the updated S 

W ^ 

values. 

7) Fresh computation is initiated for the next time step starting 
from step 2. 
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4.3 APPLICATION OF NEWTON - RAPHSON TECHNIQUE 

A Newton-Raphson technique has also been used in the present work 

to solve the system of coupled pressure equations, Equations 

4. 2-4. 5. The Newton-Raphson iteration formulae are described in 

Appendix A2. Owing to abrupt changes in the value of capillary 

pressure, P^ in the low saturation zone. Figure 4.4, calculation 

of the derivative, dS /dP , in this zone becomes highly 

w c 

susceptible to errors. Further, there are abrupt changes in the 
value of derivative itself. Hence Newton-Raphson iterations are 
not sustained and the scheme diverges in general after a few 
iterations. One has to employ very small time steps and small 
relaxation factors to sustain the iterations. Hence our study 
shows that Newton-Raphson iterations are not suitable for oil 
recovery simulation since the constitutive relations do not 
possess smooth features. However a truncated Newton-Raphson 
scheme has been successfully used in this study. Here the 
derivative terms for the rapidly changing properties are dropped 
from the elements of the Jacobian matrix. This is seen to give 
accurate results using the same computational time as in Picard 
Iterations. Since there is no additional advantage using the 
truncated Newton-Raphson method, all results have been presented 
here using Picard iterations alone. 

4.4 ENHANCED OIL RECOVERY SIMULATION USING DOMAIN DECCM^OSITION 
We now describe the oil recovery simulation using domain 
decomposition. The physical domain is shown in Figure 4. 1 and the 
initial and boundary conditions are identical to those consldred 
in Section 4.1. Uzawa’ s algorithm is used for interface 
treatment. Hence the computer code is suitable for parallel 
computing. Implicit formulation is used in each sub-domain to 
circumvent time step limitations. 
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dP/dy = 0 



Figure 4.5 

As shovm in Figure 4.5, the physical domain is split into two 
equal sub-domains. It should be noted at this stage that Uzawa’s 
algorithm can accomodate any number of sub-domains without any 
add! tonal complications. We limit ourselves to two sub-domain 
formulation to make the mathematical description more 
comprehensible. The interface, C, is a straight line parallel to 
the y-axis. This choice of interface saves a considerable amount 
of effort that may be required while calculating the derivatives 
along the outward drawn normal. Thus, for the sub-domains 1 and 
2, boundary conditions on three edges are automatically prescribed 
from the original problem and only the Interface has to be taken 
care of. We are solving for the phase pressures, and P^, in 
each sub-domain. It should be noted that the problem is 
non-linear. Hence the quantities to be matched at the interface 
are the phase pressures , P^ and P^ and the Darcy velocities of 
oil and water. This is in analogy with flux matching as described 

in chapter 3. These Darcy velocities are, Equations 2. 4-2. 5, 

K K ap K K ap 

- LI - and — Since K and p are constants with 

p ax p ax 

respect to pressure for the Newtonian case, we are left with the 
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pressure fluxes, - K 


and K 


to match. 


Derivatives 


with respect to x are calculated along the outward normal of the 
sub-domains. We now outline the algorithm developed here: 

1) At the beginning of each time step we supply the guess values 

of the pressure fluxes at the interface. The most judicious 

choice would be the fluxes calculated in the previous time step. 

Let us call them and . 

1 ? *' 

2) With X^ and X chosen, we solve the pressure equations in 
sub-domains 1 and 2 using X^ and X^ as the boimdary conditions at 
the interface. A central difference scheme as described in 
Section 4.1 is used to solve the pressure equations in regions 1 
and 2. Let the solutions be and in region 1 and P^ and 

j ol wl ® o2 

P in region 2. 

m2 


3) Using Uzawa’s algorithm we update the fl\ixes as 

X^ = X^ + p f P^ I - P^ I 1 

“ ° I oilj 

X^ = X^ + p f P^ I - P^ I ] 

W w ^ w2 1 I J 


(4.6) 


(4.7) 


A value of p = 0.25 has been found to be optimum in this study. 

4) Steps 2-3 are repeated until the conditions 


< e and 


X^^^- X^ 

w w 


(4.8) 


are satisfied. Here I Indicates the level of Interface iteration. 

5) Using the pressure values calculated in regions 1 and 2, S^ is 
updated and ppv is calculated. 

6) We proceed to the next time step and the sequence of 
operations, steps 1-5, is repeated. 

In Table 4.1 we compare the saturation values calculated 
using the full domain simulation and domain decomposition codes. 
The values are seen to be extremely close . Table 4.2 shows the 
comparison between the CPU times taken for three hour simulation 
using full domain simulation and domain decomposition code on a 
refined grid of Ax == 0.0167m. Domain decomposition can be seen to 
be faster even on a sequential machine. In Table 4.3 we show the 
effect of number of sub-domains on CPU time through a non-linear 
heat conduction problem using Gauss-Seidel iterations. 
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Table 4.1: Comparison of full domain simulation and 

simulation with domain decomposition for the oil recovery 

problem. 

X in meters and t in hours. Table shows water saturation values 
at a formation temperature of SO^C. Newtonian behaviour of oil is 
considered. 

A: Full domain simulation B: Domain Decomposition 


X =► 

0.5 

1 

1.5 

2.0 

2.5 

3.0 

3.5 

4.0 

A 

t- 1 S 

.3295 

.2353 

.2232 

.2011 

.2000 

.2000 

.2000 

.2000 

B 









.3297 

.2349 

.2232 

.2012 

.2000 

.2000 

.2000 

.2000 

A 

t= 3.0 

.3656 

.2760 

.2532 

.2325 

.2198 

.2055 

.2000 

.2000 









B 

.3656 

.2762 

.2539 

.2340 

.2202 

.2059 

.2000 

.2000 


Table 4.2: Comparison of CPU times for full domain simulation and 
simulation with domain decomposition ( two sub-domains) for three 
hour simulation of the oil recovery problem in the presence of 

surfactant. 

CPU time in seconds. Temperature T in “ C. 

Oil is taken to be Newtonian. 

Domain 

Full domain Decomposition 

T=30 
T=40 


900 

845 

920 

866 

933 

917 

— — 


T=60 










Table 4.3: Comparison showing the effect of number of sub-domains on 

CPU time. 


CPU time (in seconds) are for steady state simulation 
non-linear heat conduction problem described in chapter 
Gauss-Seldel iterations. N is the number of sub-domains. 
Indicates full domain simulation. 

A: 241 nodes 
B: 3001 nodes 



N=:0 

M = 2 

K = 3 

A 

1.0 

1.1 

0.85 

B 

23.2 

20.6 

18.4 


of the 
5 using 
N = 0 






BC ot 



Figure 4.3 Matrix structure In pressure 







L^rtAT 1 im b : KfcbULTS 


In the present chapter the results have been given under the 
following heads: 

1) Effect of surfactants on oil recovery. 

2) Effect of shear-thinning behaviour on oil recovery. 

3) Effect of enhanced inflow plane pressure on the saturation 
profile and oil recovery. 

5.1 DATA 

Fluid properies such as oil and water and formation properties are 
given in Table 5.1. The constitutive relations namely^ k = k (S ) 

r r w 

and P^= P^(S^) are shovm in Figures 5.1a and 5.1b. Oil is 
modelled as a Newtonian fluid (n=l) as well as a shear thinning 
fluid (n<l). For the study of shear thinning behaviour, modified 
Darcy’s law, Equation 2.18, is used and the value of n is taken to 
be 0.6. This value of n encompasses a variety of crude oils (.Aziz 
(1986)) that are generally found in oil reservoirs. The formation 
temperature and injection temperature of water are taken to be 
equal. Temperature is however used as a model parameter as oil 
and water viscosities are strong functions of temperature. The 
ppv values and saturation profiles presented in this chapter 
correspond to a three hour simulation. 

5.2 GRID 

The grid sizes used are Ax = 0.05m and Ay » 0.2m and the time step 
used varies from At=3.6 seconds for the initial time levels to 
At=36 seconds for the later ones. Further grid refinement does 
not show any significant changes in the numerical results and this 
point is demonstrated in Table 5.2. A convergence criterion of 

100 X ? ^ s 0.01 (5.1) 

where 0 is P or P , is used for the Picard iterations within a 

0 W 

time step. 

In order to minimize spurious oscillations, smoothened 
initial conditions for P , P and S are used. These spurious 

o w w 

oscillations are observed for short times when initial conditions 
are not smooth. Initial conditions are of the form 
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For 0 s X s 0. 4, 


P = 1.79 - 1.225X 

o 

S = 0.86 - 1.65X 

w 

P = P (S ) 

c c w 
p =p - p 

wo c 

For X > 0.4, 

P = 1.31 MPa (5.6) 

O 

= 0.2 (5.7) 

P = P (S ) MPa (5.8) 

c c w 

P = P - P MPa (5.9) 

woe 

The computer code deveoped here is written in FORTRAN 77 and has 
been implemented on HP 9000. The full domain simulation code has 
around 600 lines. The domain decomposition code has been written 
in such a manner that it can be implemented on a parallel computer 
without major modifications. Calculations in sub-domains are 
performed in separate subroutines and a main program similar to 
the front-end-processing part described in Chapter 4 invokes them 
and performs the necessary operations. The main program has been 
written in a versatile manner so that it can accomodate as many 
sub-domains as necessary. 

5.3 TESTING 

Veracity of all the computer codes develoi>ed namely, full domain 
simulation with PJeard iterations and Newton-Rapson iterations, 
domain decompostion code and the matrix inverter has been tested 
by solving the non-linear heat conduction problem 

7. K(T) V T = -U- (5.6) 

With the same Initial and boundary conditions as described in 
Appendix Al. The reference solution is taken from an independent 
computer program that uses Gauss-Seidel Iterations for the same 
problem. Values for comparison are shown in Table 5.3. The 
analytical solution for T (one dimensional problem) at steady 
state for K(T) = 1 - 0.4T is 

5 - ■/l6x-t-9 ^5 7 ) 

2 


MPa 

(5.2) 


(5.3) 

MPa 

(5.4) 

MPa 

(5.5) 
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5.4 EFFECT OF SURFACTANT ON OIL RECOVERY 

Results obtained by modeling oil as a Newtonian fluid are 
considered first. Figure 5.2 shows oil recovery in terms of ppv 
as a function of the formation temperature. The three cases 
namely, injection without surfactant, with surfactant and the 
ideal case are compared here. It is seen in Figure 5.2 that ppv 
increases with temperature. An increase in temperature lowers oil 
viscosity and hence the mobility of the oil-water bank. The 
effect of lowering surface tension in the presence of surfactants 
can also be seen in the figure. A dramatic increase in ppv of 
about seven times is seen in the presence of ideal surfactants. 
Even in the presence of partially effective surfactants the 
increase is appreciable. 

Figure 5.3, compares water saturation profiles after three 
hours of injection with and without surfactants. It can be seen 
that water migration is diffuse in the case of injection without 
surfactants and with partially effective surfactants. The 
saturation curve for the latter case is always higher than for the 
former case showing greater mobility of the oil-water bank in the 
presence of surfactants. In the Ideal case water moves as a front 
since capillary resistance is zero. Surfactants that show this 
degree of effectiveness remain to be synthesized at this time. 

5.5 EFFECT OF SHEAR THINNING BEHAVIOUR OF OIL 

Figure 5.4, delineates the effect of shear-thinning behaviour of 
oil. Oil recovery in the presence of shear-thinning behaviour is 
compared with that in the presence of Newtonian behaviour for all 
the three cases mentioned before. Shear-thinning behaviour is 
characterized by a decrease in the slope of shear stress-shear 
strain curve at increasing values of shear stresses. Accordingly 
an increase in oil recovery is to be expected. This is seen in 
Figure 5.4. 

5.6 EFFECT OF INFLOW PLANE PRESSURE 

The effect of enhanced inflow plane pressure is presented in this 
section. A higher pressure value of 3.63 MPa Instead of 1.79 Mpa, 
that is used throughout this analysis, is used on the inflow 
plane. Formation temperature used is 30 C and fluid injection is 
in the presence of surfactant. Newtonian behaviour of oil is 
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consiared here. With a 4m long domain the saturation profile 
shows non-physical oscillations and after 20 minutes it breaks 
down. This may be attributed to the proximity of the outflow plane 
that precludes the use of a Dirichlet boiandary condition. As 
discussed in Chapter 2, the use of Dirichlet boundary condition on 
the outflow plane is justified only if the saturation profile in 
the vicinity of the outflow plane is not greatly disturbed. When 
the inflow plane pressure is raised substantially the water-oil 
interface moves at a greater speed and reaches the outflow plane 
well before three hours. Therefore the saturation profile breaks 
down soon since the use of Dirichlet boundary condition is not 
proper. An interface of two immiscible fluids like oil and water 
is inherently unstable. Therefore with increasing inflow plane 
pressure the interface, which as such is unstable, becomes more 
prone to a collapse. 

With the same enhanced inflow plane pressure, with a 8m 
long domain the saturation profile is continuous even after three 
hours. This simulation uses the domain decomposition code with 
two 4m long sub-domains. Results are shown in Table 5.4. 
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Table 5.1 ; Properties Used in Numerical Calculation 


Absolute Permeability, K = 132 Darcies 

Porosity, e = 0.375 

Compressibility: oil, ^ = 0,03447 Pa"^, water, f = 0.02137 Pa 

o 

In-situ formation pressure = 1.31 MPa 

Injection pressure = 1.793 MPa 


Dynamic Viscosity (Pa-s) 




T, °C 20 

40 

60 

80 


Oil 

4.12 

2.37 

0.61 

0. 10 

Water 

1.03xl0~^ 

6.6x10“'^ 

4.73x10"'^ 

3.67x1 



J.OU 11 snowing the effect of grid refinement on 

saturation values in the oil recovery problem (Newtonian behaviour 
of oil) in the presence of surfactant.. Formation temperature is 
40°C. 

A: Ax = 0.05m ; B: Ax = 0.033m 
Comparison of saturation values. 


X ^ 

0.5 

1 

1.5 

2.0 

2.5 

3.0 

3.5 

4.0 

A 

t= 1 0 

.3164 

.2435 

.2156 

.2040 

.2000 

.2000 

.2000 

.2000 

B 









.3146 

.2461 

.2084 

.2051 

.2000 

.2000 

.2000 

.2000 

A 

t» 3.0 

.4015 

.3376 

.2786 

.2520 

.2212 

.2110 

.2000 

.2000 









B 

.4021 

.3382 

.2772 

.2516 

.2219 

.2106 

.2000 

.2000 


Comparison of ppv values. 


1 

2 

3 

4.0616 

7.8607 

11.6167 

4.0594 

7.7132 

11.4211 


B 











Table 5.3: Comparison of the solutions of the non-linear heat 

conduction problem ( chapter 5) obtained by using direct matrix 
inverter and Gauss-Seidel iterations. 

A: direct matrix inverter 
B: Gauss-Seidel Iterations. 

All quantities are dimensionless. 


X ^ 

0.2 

0.4 

0.6 

0.8 

1.0 

2.0 

4.0 

A 

4. — . IQ 

.6634 

.3872 

.1983 

.0893 

.0356 

.0001 

0.0 

* lo 

B 

.6633 

.3872 

.1974 

.0864 

.0355 

.0001 

0.0 

A 

j. — n c: 

.8146 

.6395 

.4830 

.3504 

.2439 

.0209 

0.0 

B 

.8146 

.6400 

.4831 

.3504 

.2438 

.0209 

0.0 


Table 5.4: Results of enhanced Inflow plane pressure simulation 
of 8m long domain with two 4m long sub-domains. Duration of the 
simulation is three hours and fluid injection is in the presence 
of surfactant. Oil is taken to be Newtonian. Formation 
temperature is 30°C. 

Inflow plane pressure: 3.63 MPa 
X in meters and t in hours. 

X 

S 


1 2345678 

0.345 0.272 0.265 0.241 0.212 0.202 0.200 0.200 


t 

ppv 


0.5 


1.0 1.5 2.0 2.5 3 


2.217 3.770 5.310 7.001 8.752 10.422 











'function of formation temperature 







Figure 5.4 Effect of shear-thinning behaviour of oil. 
surfactant b) With surfactant c) Ideal case 


a) Without 



CHAPTER 6: SCOPE FOR FUTURE WORK 


The following extension to the present work is recommended to make 
numerical simulation of oil recovery complete and realistic. 

1) Radial simulation: 

In this work we have considereda laboratory scale rectangular 
domain with one injection well and one production well. In 
reality, reservoirs are characterized by one centrally located 
injection well radially sorrounded by production wells. Analysis 
of these reservoirs involve the use of radial coordinates. 
Instead of performing the analysis for the whole system, which may 
be very time consuming, a novel approach would be is to assume 
symmetry and analyze only one pair of injection well and 
production well and extend the results to the whole domain. In 
the vicinity of the wells near-field effects are predominant and 
they should be tackled carefully. Domain decomposition finds 
utility here as the affected regions can be localized by assigning 
sepa,rate sub-domains to them and can be studied closely. 

2) Stabllty limits of simulation: 

The numerical simulation presented in this work has been done for 
one particular formation l.e. one particular set of formation and 
surfactant properties. Owing to high degree of non-linearity 
involved in the governing reservoir equations, applicability of 
this simulation code to other sets of properties cannot be 
guaranteed. Hence to ascertain versatility of the present 
formulation and computer code, simulation should be performed with 
other formation properetles as well. There is considerable 
evidence that with certain properties the water saturation 
profiles become unstable. This leads to mixing of oil and water 
and a fall in oil recovery with increasing injection pressure. 
This is a practical problem that deserves attention. 

3) Three-phase flow: 

Besides oil and water, the third phase that is present in certain 
reservoirs is air, though comparatively much less in quantity. A 
three phase simulation (oil, water and air) would thus lend more 
credibility to this work. An oil-water-air-steam simulation. 
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t-hougii dittlcult inust also be taken up for analysis. 

4) Adsorption-Desorption In the formation; 

In the case of surfactant flooding a considerable amount of 


surfactant 

is lost to the formation 

because 

of adsorption. 

Desorption 

i.e. flow of 

surfactant 

from formation to 

the 

surfactant 

solution also 

takes place, 

though 

much less 

in 

magnitude. 

Bacause of 

these mass 

transfer 

phenomena 

the 


constututlve relations must be modified in conjunction with the 
variation of concentration of surfactant solution. An additional 
mass transfer equation comes into picture and it must be solved 
along with the pressure equations. 

5) Three dimensional simulation; 

An extension of the present is a three dimensional simulation 
including Inhomogeneity and anisotropy of the formation. This 
will make reservoir simulation computationally more expensive. 
Domain decomposition with parallelization developed in the present 
work is the best option to tackle this problem. 

6) Parallelization; 

Implementation of the parallel code developed here on a parallel 
computer remains incomplete. It could not be done because of lack 
of availlbility of parallel computing facility at IIT Kanpur. 
This implementation would not only save a lot of computational 
time but it would also enable us to perform simulation for real 
life reservoirs. 

7) Interface treatment; 

Uzawa's algorithm that is used here is one of the rudimentary 
techniques for interface treatment in domain deconposltion. Many 
advanced parallelizable techniques are now available (GJowinski 
(1983)). Use of these techniques can be expected to enhance the 
convergence rate of the Interface iterations and make the whole 
simulation faster. These techniques can be tested systematically 
in the future. 
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APPENDIX Al: ITERATIVE TECHNIQUES FOR NON-LINEAR EQUATIONS 
In this section we describe two iterative schemes, Picard 
iterations and Newton-Raphson iterations, that have been used in 
this study to solve non-linear algebraic equations. 


I 


Al.l PICARD ITERATIONS 

Let us first consider the single variable problem. The equation 
to be solved is 

F(x) * 0 (1) 

where F(x) is a non-linear function of the variable x. In this 

technique, we first recast Equation 1 in the following form 

X = f(x) (2) 

This recasting can be done in many ways and there is no uniqueness 

in this step. One then assumes a starting (or first iterate) 

value, x^, for the variable x and uses this in f(x) of Equation 2 

2 

to get a new and hopefully, improved estimate, x . Thus 

x^« f(x*) (3) 

x^ is now used as the argument of f(x) in Equation 2 to yield the 
next estimate, x^. This procedure is continued and we have the 
algorithm as 

x''*^= f(x’') ; k - 1, 2, ... C4) 

Computations are stopped when the successive Iterates are such 


that 



a e 


(5a) 


and/or 

\F{x'^*^)\sc (5b) 

where e is the desired tolerance. 

This technique can now be easily extended to the more general 
case Involving N variables, x^, x^, ..., x^. The i equation , 
F (x) =0, in Equation 1 can be rearranged as 

x^ = f^(x) ; 1 “ 1, 2, ...» N (5) 

In the matrix form it may be written as 



"'l 1 


f^(x) n 

X - 


= f(x) = 

f^(x) 


- - 




(7) 
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In analogy with Equation 4, the algorithm can be writen as 

k 




k , 


( 8 ) 


One starts with a choice , x^, of N variables and iterates until 
convergence criteria similar to Equations 5a and 5b are satisfied 
on all x^ and F^[x]. 


A1.2 NEWTON - RAPHSON ITERATIONS 

A very popular method for solving sets of non-linear algebraic 
equations is the Newton-Raphson technique. 

For the single variable case, we have 

F(x) = 0 (10) 

As described before, x*^ is the estimate of the root, a, at the 

k^l 

iteration and x is the improved estimate at the next iteration, 
then the Newton-Raphson algorithm is 


k+l 

X 


k 


X 


FCx*^) 
f’ (x*') 


( 11 ) 


* k dF* k 

where F (x ) ( = ^ ) is evaluated at x . The convergence 

criterion of the Newton-Raphson technique is that if the function, 

F(x), has a continuous derivative in the neighborhood of a root, 

> 

a, and if the derivative, F (x), is non-zero in the neighbohood of 
this root, then the Newton-Raphson technique converges to a, 
provided x^ (the first guess) is taken sufficiently close this 
root. Computations are stopped when the conditions in Equations 
5a and 5b are met. 

The extension of this technique to a set of N equations in N 
variables, x , x , . . . , x , is straightforward. The i^^ equation 

X 2 N 

from the set. Equation 10, is 

F (x) = 0 = F (x , X , ..., x„) ; 1 = 1, 2 N (12) 

1 1 X 4 Q| fl 

Iterations are then governed by 


F^[x‘] . 
F^[x‘'] 


aF^/3X^ 

SF /ax 

2 1 

aF /ax, . 
1 2 

aF/ax^ . 


X - X ^ 

^2 - ^ 

F [x*'] 

W 

NXl 

aF /ax 

■- N 1 

aF /ax . 

H 2 

, aF /ax 

N N 

KXN 

* k 

Lx «. X -* 

U 


( 13 ) 
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or 


- Ft x‘'] = At t x^*^- x*'] (14) 

where At x*'] is the Jacobian matrix evaluated at x*' and it is a 
matrix of constants. Equation 14 can be rewritten as 

= x*' - 1^ At x'' ] j Vt x’' ] (15) 

It can be easily seen that Equation 11 Is a special case of the 
more general Equation 15. Iterations stop when the convergence 
criteria as outlined in Equations 5a and 5b are satisfied for all 
Xj and Fj tx] . 


A1.3 APPLICATION OF NEWTON - RAPHSON TECHNIQUE 

We now use the Newton-Raphson technique to solve the 

dimensional non-linear transient heat conducton problem. 

„ „ „ „ dT 

X.' ^ I ^ ^ -w~ 


two 


(16) 


with the 

initial 

and boundary conditions. 



at t=0 

T = T 

o 

at 

t>0, 

at X = 0 

T «= T 

1 



at X = X 

T = ■ 



at y = 0 

3T 

dy 



at y = Y 

dT 

dy 


(17) 

(18) 

(19) 

( 20 ) 

( 21 ) 


In Equation 16, K is the thermal conductivity and its variation 
with temperature is given as 

K = 1 - pT (22) 

Equations 16-21 are discretized using central differencing with 
implicit time stepping. The control volume is shown in Figure 
Al.l. 



43 



t and J are the grid locations in x and y directions respectively. 
N is the number of nodes in the x-direction and M that in the 
y-direction. Ax and Ay define a Cartesian mesh and At is the time 
step. (n+1) indicates the current time step and n the previous 
time step. Upon discretization, we get 
For the interior points; i = 2, N-1 and J =2, M-1 




and the discretized boundary conditions are: 


for i = 1, 
for i = N, 
for J = 1. T 
for J = M. T 


T * T 

i.J 1 

T * T 

i.J 2 

- T = 0 
i,J+i i.J 

- T = 0 
i.J i,J-l 


(24) 

(25) 

(26) 
(27) 


K , K , K , K are the harmonic averages of the values of K at 

e K n 8 

nodes that lie on either side. Thus, by using finite differences 
we have reduced the governing partial differential equation and 
the boundary conditions. Equations 16-21, to a system of NM 
non-linear algebraic equations, Equations 23-27, which may be 
written as 




12 ’ 


T , T ] 

H,M-1 MM 


(28) 


where the NM unJfcnowns are the temperatures at the NM grid points. 
For simplicity, we reduce the two dependent subscripts, i and J, 
to one subscript, J + (t-l)M, where i and J hold the same meaning. 
Equation 28 may now be written as 

(29) 

for all i and 


F [ T ] = 0 


where the NM unknown temperatures are now T 

J. 




The system of equations. Equations 23-27, is now solved using 
the Newton-Raphson technique. The first step would be is to form 
the Jacobian matrix A. A close look at Equations 23-27 shows that 
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any F has only the temperature at the grid point i.j and the 
temperatures at the four neighbouring grid points as its 
variables. Therefore, we evaluate only the derivatives with 
respect to these temperatures and safely put other derivatives to 
zero. The non-zero derivatives are 

3F , 

J+(l-l)M 

1) gj 

T 

I.J 

For i = 2, N-1 and J » 2, M-1 


r ^ 

k k k 


T r 

, ak 

1 • 

+ ** + “ + " 

+ ^ 

1 T f 

1 e 

L Ax^ 

.2 .2 .2 

Ax Ay Ay 

At 


A 2 aT 

Ax 1 , J 

ak 

1 1 
+ L + 1 

ak T 

T 

ak T 


■ 1 

^ 1+1, J 

e 1-: 

aT 

. 2 0T . 2 

aT 

. 2 

aT . : 

I.J 

Ay 1 , J Ay 

I.J ^ 

Ax 

I.J Ax 


ak 


T ak 

1,J+1 n 


T ak 

l.J-l 8 


. 2 aT . 2 aT 

Ay i.j Ay 1 ,J 


For i = 1 =1 

For i = N =1 

For J = 1 = -1 

For J = M =1 


2 ) 


aF 

aT 


T 


1+1, J 


(31) 


For i * 2, 



For i = 1 
For i = N 
For j = 1 
For J = M 


1 + 1 , J Ax' 

= 0 
= 0 
= 0 
= 0 


N-1 and j * 2, 



M-1 

T 

1+1. J 


ak 


1+1. J 


(30) 
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ar 


3) 


J+ (l-l) M 


ai 


1-1. j 

For i = 2, N-1 and J 


2, M-1 


T ak 

1, J w 


Ax 


ax 


+ — + 
2 


T ak 
1-1. J *• 


1-1, J Ax 


Ax 


ai. 


1-1, J 


For 1 = 1 
For i = N 
For J = 1 
For J = M 
aF 


= 0 
= 0 
= 0 
= 0 


4) 




ax 


i.j+i 


For 1=2, N-1 and J = 2, M-1 


ak 

n 


+ -ii- + 


I, J 

A 2 5X .2 

Ay 1,J4-1 Ay 


X ak 

i.J*i n 


A 2 ax 

Ay i,j*i 


For 1 
For i 
For j 
For J 


1 

N 

1 

M 


= 0 
= 0 
= 1 
= 0 


aF 


5) 


J+(1-1)H 


ax 


i, j-i 


For 

i = 

2, N-1 

and j = 


X 

ak 

k 

_ _ 

1, J 

s 

_ + _!L 


. 2 

ax 

. 2 


Ay 

1, J- 

-1 Ay 

For 

i = 

1 

= 0 

For 

i = 

N 

= 0 

For 

J = 

1 

= 1 

For 

J = 

M 

= 0 


X ak 

1. j-i » 

. 2 ax 


1, j-i 


(32) 


(33) 


(34) 
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In all the derivatives calculated above, derivatives of 
Kj(i=e,w,n,s) with respect to T are calculated along the same 
lines as the example given below. 


dK 


dK 


dX 


= and 


dl 


= 0 


(35) 


1 + 1. J “ 1-1, J 

With the derivatives calculated, the Jacobian matrix A is formed 
as described in Equation 13. We now outline the algorithmic steps 
to solve this transient problem. 

1) For every time step choose an initial guess (first iterate) for 
the temperatures. This initial guess may be readily chosen as the 
temperature values at the previous time step. 

2) With the initial guess chosen, iterate the temperatures within 
the time step using the Newton-Raphson scheme, Equation 15, 


[ I = [ T I*' - A [ 1“] j F [ t‘'] (36) 

where A [I*'] is the Jacobian calculated at [ T Similarly 

F [ T*'] are the functions, F calculated at [ T ]*'. 

3) Continue the iterations until convergence is achieved. The 
convergence criterion can be written as, 

I [ T - t T ]“! s e (37) 

where e is the desired tolerence. 

4) Upon convergence proceed to next time step and go through the 
same sequence. 
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APPENDIX A2: MATRIX STRUCTURE FOR OIL RECOVERY PROBLEM 


In this section we list the values of coefficients of the 
discretized pressure equations, Equations 4.2 and 4.3. These 
coefficients are the non-zero entries of the matrix described in 
Chapter 4 and shown in Figure 4.3. 

A2.1 PICARD ITERATIONS 


On using Picard iterations, the coefficients are as follows: 
k k_ p k k p 


Here 0 = ° and y = 


G u 


e |i_ 


p (Ax)^ /i,j 


t = [■ 


dS 


r 

Hh 

1 

/ 0 + 0 

1 * " t 

0 + 0 \ *1 
n si 

L At 


^0 

^ (Ax)^ 

(Ay)^ A, j J 


n+1 


f dS 

i “ 

s n+ 1 

1 


1 

II 

O 

Q 

/ 0 

f « 1 

1 dP 

' c 

M. J 

At 

» 

^ p (Ax)^/ 1 , j 
0 


/ 0 \ n4- 1 

^ p (Ay)^ / 1 , j 

O 


dS 


S 


dP 


G = 


V V** 

r / 

\ n-f 1 

P'‘ + i- I 

dS \ n+l 1 

"1 P" 

Li At 

/ 

• o At \ 

dp, ; « J 






p (Ax) 

w 

dS 
S C 

vr M dr 


B 


At 


_1_ 

P 


I 


r + V 

e H 


7 + r 


1 


i (Ax! (Ay^ 


n+l 
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In expressions for the coefficients A,B...G given above, the 
suffixes e, w, n and s stand for the locations midway between 
(i,i+l), (i,i-l), (j,J+l) and respectively. The 
corresponding properties are calculated as the harmonic average of 
their values prevailing at the nearest nodes. For example, = 
harmonic average of (yCi+l), y(i)). 


A2.2 NEWTON - RAPHSON ITERATIONS 


For Newton-Raphson iterations the matrix is the same as in Picard 
iterations. The non-zero entries i.e. A-G, A-G are as given 

O 0 w w 

below. Since P = P - P , it should be noted that 

COW 


ao 

ap 


ao 

dP 


aQ_ 

ap 


where Q may be any quantity. 


Further, the derivatives of 0^ (i = e, w, n and s) and (i = e, 
w, n and s) with respect to P and P are calculated on the same 

O W 

lines as the following example; 


a 0 2 k 

O ^ 


r 1 ^ j 

r' \ 

1 1 

ap “ SF 

oH-l, J 

*^0 

k k 

L j rol , J 

J i 

k 

rol+l,J J 


, dk ^ 


It should be noted that 0^ is the harmonic average of e. . and 0;^, 
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other derivatives of 0 and y are calculated similarly. Here 


%= ( ■ V ( 

C j 


~ ■ dP 




dS 

and b = b = 

o w dP 

c i, j 
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A = - -ii 

O 


r 30 

0 

p , . r 

30 T 


+ * 

ol + l.jf 
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3P " At At 3P J .2 5P , . 

wl , J J L wl , J Ax p wi , j 


p 3a 

o o 
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